# -*- coding: utf-8 -*-
"""
Created on Tue Apr  2 16:07:55 2019
周志华《机器学习》第3章 习题3.3 编程实现对率回归，并对3.0α数据进行计算
@author: songlin liang
"""

from numpy import *
import matplotlib.pyplot as plt

'''
=================================
         定义对率回归函数
=================================
'''
def LogisticRegression(X,Y,coe=1E-5):
    # 牛顿法对率回归
    # 输入X和Y，维度分别为m*n和m*1
    # coe为终止条件，当梯度小于coe时停止迭代
    # 返回拟合到的参数theta,其中最后一个参数theta[-1]为b
    m,n=shape(X)
    X1=c_[X,ones(m)]    #用c_函数实现数组的横向拼接
    theta=zeros(n+1)    #在matlab中这样生成的是 (n+1)*(n+1)的矩阵，
                        #但是在numpy中为1*(n+1)的向量
                        #注意这里向量不分行向量和列向量的
    dL=ones(n+1)        #梯度，初始任意设置为较大向量
    while linalg.norm(dL)>coe:
        p1=1/(1+exp(-dot(X1,theta)))          #计算p1,对应于(3.23)式，也就是logistic函数
        dL=-dot(X1.T,Y-p1)                    #梯度,结果为(n+1)维向量
        H=dot(dot(X1.T,diag(p1*(1-p1))),X1)   #hessian矩阵，结果为(n+1)*(n+1)
        theta-=dot(linalg.inv(H),dL)          #更新theta
    return theta                              #返回最佳参数

'''
=================================
         主程序部分
=================================
'''

#------ 西瓜3.0α 样本数据----------
X=array([[0.697,0.46],[0.774,0.376],[0.634,0.264],[0.608,0.318],[0.556,0.215],
   [0.403,0.237],[0.481,0.149],[0.437,0.211],[0.666,0.091],[0.243,0.267],
   [0.245,0.057],[0.343,0.099],[0.639,0.161],[0.657,0.198],[0.36,0.37],
   [0.593,0.042],[0.719,0.103]])
Y=array([1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0])

#--------样本数据可视化------------
plt.xlabel('x1')
plt.ylabel('x2')
plt.title('LogisticRegression on 3.0α')
plt.scatter(X[Y==1,0],X[Y==1,1],color='k',s=90)  #画散点用.scatter
plt.scatter(X[Y==0,0],X[Y==0,1],color='r',s=90)
plt.legend(['Good','Bad'])

#-------对率回归，并作图------------
# 在同一个py中调用函数LogisticRegression，需要将函数定义写在前面
#>>>>>>>>>对率回归并获取参数
theta=LogisticRegression(X,Y)
t1,t2,b=theta
#>>>>>>>>>画直线
x1=min(X[:,0])    #直线左端点横坐标
x2=max(X[:,0])    #直线右端点横坐标
y1=-(t1*x1+b)/t2  #直线左端点纵坐标
y2=-(t1*x2+b)/t2  #直线右端点纵坐标
plt.plot([x1,x2],[y1,y2])         #画直线用.plot
#>>>>>>>>>画正负号标记
xm=(x1+x2)/2
y_positive=(1-(t1*xm+b))/t2
y_negative=(-1-(t1*xm+b))/t2
plt.text(xm,y_positive,'+',fontsize=30)
plt.text(xm,y_negative,'-',fontsize=30)
#>>>>>>>>>写入对率回归得到的参数
plt.text(x1,y1,('%.2f x1%+.2fx2%+.2f=0'%(t1,t2,b)))